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ABSTRACT 

We  present  the  numerical  solution  to  a  problem  of  maximizing  the  lift  to  drag  ratio  by 
rotating  a  circular  cylinder  in  a  two-dimensional  viscous  incompressible  flow.  This  problem 
is  viewed  as  a  test  case  for  the  newly  developing  theoretical  and  computational  methods 
for  control  of  fluid  dynamic  systems.  We  show  that  the  time  averaged  lift  to  drag  ratio 
for  a  fixed  finite-time  interval  achieves  its  maximum  value  at  an  optimal  rotation  rate  that 
depends  oi\  the  time  interval. 
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1.  INTRODUCTION 


Active  control  of  fluid  dynamic  systems  is  an  area  of  research  that  has  been  con¬ 
siderable  growth  during  the  past  ten  years.  This  activity  is  motivated  by  several  potential 
applications  in  engineering  sciences  and  optimal  design.  The  development  of  theoretical 
frameworks  and  the  construction  of  practical  computational  algorithms  for  control  design 
are  highly  complex  problems.  In  this  paper  we  consider  the  problem  of  maximizing  lift  to 
drag  ratio  by  rotating  a  circular  cylinder.  This  problem  is  simple  enough  to  be  solved  by 
direct  numerical  calculation  (as  we  do  here)  and  complex  enough  to  serve  as  an  excellent 
test  problem  for  future  theoretical  and  computational  approaches.  The  main  objective  of 
this  short  note  is  to  provide  the  solution  to  two  optimization  problems  that  can  be  used  as 
test  examples  by  researchers  in  this  area. 

2.  MATHEMATICAL  FORMULATION  AND  NUMERICAL  METHOD 

In  this  paper  we  consider  an  optimal  control  problem  for  a  two-dimensional  viscous 
incompressible  flow  generated  by  a  circular  cylinder  started  impulsively  into  a  combined 
steady  rotatory  and  rectilinear  motion.  This  problem  is  investigated  numerically  by  solv¬ 
ing  a  velocity/ vorticity  formulation  of  the  Navier-Stokes  equations  [2,  4].  A  nonrotating 
reference  frame  translating  with  the  cylinder  is  employed,  and  the  cylinder  rotates  in  the 
counterclockwise  direction  with  angular  velocity  Q.  The  Reynolds  number  Re  =  2 Ua/u  is 
based  on  the  cylinder  diameter  2 a  and  the  magnitude  U  of  the  rectilinear  velocity.  The 
angular /rectilinear  speed  ratio  a  =  fla/U  is  the  primary  control  parameter  in  this  paper. 
The  speed  ratio  a  can  impose  significant  influences  on  the  characteristics  of  resulting  flow 
field  as  well  as  the  temporal  evolution  of  forces  on  the  cylinder  surface. 

In  this  formulation,  the  dimensionless  governing  equations  consist  of  the  vorticity  trans¬ 
port  equation 

t  +  S-V"=iV2“  (!) 
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and  the  vector  Poisson  equation  for  the  velocity 


V2u  =  - V  x  (we*).  (2) 

The  cylinder  radius  a  is  used  as  the  length  scale  while  a/U  is  used  as  the  time  scale.  The 
dimensionless  boundary  conditions  for  the  problem  of  a  rotating  cylinder  are 


{■u  =  —a 
u  =  ex 


=  —a  sin  6tx  +  a  cos  0ey  on  the  cylinder  surface 
=  ex  at  infinity. 


This  velocity/vorticity  formulation  is  especially  well  suited  to  treating  initial  development 
of  the  flow  generated  by  impulsively  started  bodies,  in  which  the  flow  field  is  composed  of 
a  relatively  small  vortical  viscous  region  embedded  in  a  much  larger  inviscid  potential  flow. 
Consequently,  the  computational  domain  may  be  restricted  to  a  smaller  region  where  the 
vorticity  contributions  are  contained. 

The  vorticity  transport  equation  (1)  is  first  discretized  by  a  second  order  central  difference 
in  the  radial  direction  and  a  pseudospectral  transform  method  in  the  circumferential  direction 
for  all  spatial  derivatives.  This  semi-discretization  form  of  equation  (1),  consisting  of  a  system 
of  ordinary  differential  equations  in  time  can  written  as 

—  =  F(u),  do  =  (u>2i2,  •  •  •  (4) 

for  all  the  interior  grid  points.  Therefore,  the  calculation  procedure  consists  of  the  following 
steps  to  advance  the  solution  for  any  given  time  increment: 

Step  1:  Internal  vorticity  over  the  fluid  region  at  each  interior  field  point  is  calculated  by 
solving  the  discretized  vorticity  transport  equation.  An  explicit  second-order  rational  Runge- 
Kutta  matching  scheme  based  on  the  work  of  [8]  is  used  to  advance  in  time  for  (f). 

The  discretization  in  time  of  (4)  thus  can  be  written  as 

-n+l  -  n  ,  (&,&)-  Mu§\) 

W  W  +  (5) 


g\  =  F(«")Ai, 
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92  =  F(un  +  0.5g\)At, 

a  n  A  A 

93  —  *9 1  —  9 2> 


(6) 


where  (<fi,  #3)  denotes  the  scalar  product  of  <fi  and  53.  This  step  consists  of  the  kinetic  part 
of  the  computational  loop. 

Step  2:  Using  known  internal  vorticity  values  at  all  the  interior  grid  points  from  step  1 , 
the  generalized  Biot-Savart  law  of  induced  velocity 

<2(r,t )  x  (f  -  fl) 


=  -hJL^ S 


-dA 


-h!L 


2Ul{f,t)  x(r-rb) 


n 2 


<L4-f  £/ 


(7) 


is  wsed  to  update  the  boundary  vorticity  values  at  all  the  surface  nodes.  Here  fi  represent  all 
grid  points  located  on  the  solid  boundary. 

One  of  difficulties  encountered  in  the  simulation  of  viscous  flow  is  to  prescribe  the  ap¬ 
propriate  nonvelocity  boundary  conditions  at  the  solid  surface.  The  boundary  vorticity  is 
required  for  the  formulation  based  on  the  velocity/vorticity  (or  stream-function/vorticity) 
variables.  In  order  to  overcome  this  particular  difficulty,  we  pose  the  kinematic  relation¬ 
ship  between  velocity  and  vorticity  fields  on  the  infinite  domain  in  the  integral  form  of  (7). 
This  boundary  integral  method  proposed  by  Wu  and  Thompson  [9]  provides  the  basic  link 
between  the  velocity  and  vorticity  fields  throughout  the  numerical  computations. 

Step  3:  At  this  stage,  all  the  vorticity  values  in  the  computational  domain  are  known  at 
the  new  time  level.  Then,  the  velocity  at  points  on  the  outer  perimeter  of  the  computational 
domain  is  calculated  by  this  integral  kinematic  constraint 

u(r,t)  x  (r-  ro) 


^  1  f  f  W(M)  x  (r  ■ 

-  -TJ\d - 


-dA 


-rJJB 


2f i(f,  t)  x  (r  —  ro) 


'-rol2 


dA  +  U. 


(8) 


Here  fo  denotes  the  points  located  on  the  outer  perimeter  of  the  computational  domain. 

In  fact,  this  integral  representation  allows  us  to  determine  the  velocity  point-by-point 
explicitly  if  all  vorticity  values  are  known  everywhere  in  the  domain  of  interest.  Moreover, 
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it  often  exhibits  more  realistic  behavior  at  the  outer  perimeter  of  the  computational  domain 
than  asymptotic  techniques  used  in  other  formulations  [1].  This  indicates  that  the  difficulty 
resulting  from  the  imposed  far-field  condition  is  removed  by  the  application  of  this  integral 
constraint. 

Step  4:  The  new  velocity  field  can  be  established  by  solving  the  Poisson  equations  (2)  with 
prescribed  solid  boundary  conditions  and  outer  boundary  conditions  that  have  been  determined 
from  step  3. 

The  final  form  of  the  discretized  Poisson  equations  can  be  written 

{£:l  <•> 

where  u  and  v  are  the  vectors  of  unknown  interior  nodal  values.  The  resulting  11-banded 
matrix  equations  are  then  solved  by  a  preconditioned  biconjugate  gradient  routine  [3].  This 
step  completes  the  computational  loop  for  each  time  level. 

One  further  important  point  to  be  noted  in  the  boundary  integral  approach  is  the  deter¬ 
mination  of  the  initial  flow  field.  In  contrast  to  the  special  technique  used  by  Badr  &  Dennis 
[1],  this  integral  approach  enables  the  numerical  code  to  generate  the  initial  velocity  field  by 
executing  one  cycle  of  the  solution  procedure  (from  step  2  to  step  4)  rather  than  employing 
any  additional  treatments. 

3.  RESULTS  AND  DISCUSSIONS 

To  assess  the  numerical  algorithm,  calculations  were  performed  over  a  wide  range  of 
angular/rectilinear  speed  ratio  a  up  to  3.25  at  a  Reynolds  number  of  200.  In  this  model, 
the  rectilinear  velocity  is  fixed  as  a  constant  value  while  the  angular  velocity  is  treated  as  a 
control  variable.  Although  the  angular  velocity  can  be  time-dependent,  in  this  paper  angular 
velocity  has  been  restricted  to  the  constant  values.  For  the  works  of  a  rotating  cylinder  un¬ 
dergoing  various  kind  of  time-dependent  rotation  rate  (which  include  time-harmonic  rotatory 
oscillations  and  time-periodic  rotations  with  different  angular  amplitudes  and  frequencies), 
the  reader  is  referred  to  Ou  [6]  and  in  the  forthcoming  paper  [7]. 
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We  have  tested  the  model  against  the  experimental  work  of  Coutanceau  &  Menard  [5] 
with  excellent  agreement.  Selected  instantaneous  streamline  plots  for  constant  value  of  speed 
ratio  at  a  =  2.07  are  presented  in  figures  l(a,b).  We  note  that  for  higher  values  of  a  vortex 
shedding  continues  to  occur  as  a  result  of  cylinder  rotation  [4],  However,  the  formation 
of  the  vortex  street  behind  the  rotating  cylinder  at  high  a  (a  =  3.25)  is  not  observed  in 
the  experiments  of  Coutanceau  &  Menard.  The  experiments  were  unable  to  detect  this 
vortex  shedding  because  of  the  length  limitation  of  the  water  tank  and  the  flow  visualization 
techniques  used  in  their  work  [5]. 

The  lift  and  drag  coefficients  can  be  calculated  in  r-9  coordinates  as 

c‘4r[-i),+fw  (io> 

and 

CD  =  Ter[{^)„-Ut}siam  (n) 

where  the  subscript  b  denotes  quantities  evaluated  on  the  cylinder  surface.  In  particular, 
we  denote  the  positive  values  of  Cl  in  the  -y  direction.  The  effect  of  the  speed  ratio 
(0  <  a  <  3.25)  on  the  temporal  evolution  of  lift/drag  coefficient  is  shown  in  Figures  2. 
Note  that  the  lift/drag  increment  is  manifested  timewise  for  the  speed  ratio  up  to  a  =  2.07. 
However,  it  does  not  imply  that  as  the  speed  ratio  further  increases,  the  flow  field  will  result 
in  a  further  improvement  of  lift  and  reduction  of  drag  on  the  cylinder.  On  the  contrary,  the 
adverse  effect  of  speed  ratio  on  the  lift/drag  is  quite  evident  if  a  comparison  is  made  between 
a  =  3.25  and  a  =  2.07.  Notice  that  for  all  a  less  than  2,  in  the  initial  time  interval,  the 
slope  of  the  lift/drag  curve  seems  to  increase  with  increasing  of  a.  When  a  becomes  greater 
than  2,  the  slopes  decrease  gradually  as  a  increases.  For  all  speed  ratios  considered  here, 
a  significant  increase  in  ( CT/CdW*  is  obtained  with  increment  of  a.  However,  (Cl/Cd) 
achieves  its  maximum  value  at  much  later  time  for  a  higher  a.  We  are  interested  in  the 
problem  of  maximizing  the  averaged  {Cl/Cd)  on  a  fixed  finite  time  interval. 

Figure  3  shows  the  effect  of  the  speed  ratio  on  time-average  (for  0  <  t  <  24)  lift, 
drag  and  lift/drag  coefficients  in  the  range  of  0  <  a  <  3.25.  It  illustrates  that  the  lift 
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average  is  almost  linearly  proportional  to  the  speed  ratio.  On  the  other  hand,  the  drag 
average  remains  almost  constant  up  to  a  =  2,  then  monotonically  increases  with  speed  ratio 
thereafter.  Consequently,  the  resulting  lift/drag  average  is  not  linearly  proportional  to  the 
speed  ratio  a.  Moreover,  the  maximum  lift/drag  average  occurs  approximately  at  a  =  2.38, 
and  it  represents  a  substantial  increase  of  440%  with  respect  to  a  =  0.5. 

It  is  also  of  interest  to  study  effect  of  speed  ratio  on  the  variation  of  the  (total  lift) / (total 
drag)  force  ratio.  As  shown  in  Figure  4  this  ratio  is  similar  to  the  lift/drag  average  in  Figure 
3  and  achieves  its  maximum  value  between  a  =  2.0  and  a  =  2.38. 

4.  CONCLUSIONS 


We  have  presented  numerical  solutions  of  two  ’’flow  control  problems”.  The  first 
problem  is  to  find  a*  that  maximizes  the  time-averaged  lift  to  drag  functional 


*«*>  -  hr  [n^]  * 


and  the  second  problem  is  to  find  a £  that  maximizes  the  total  lift  to  total  drag  functional 

SlCj^dt 

M  !?CD(t,a)dt • 

Based  on  a  velocity/vorticity  formulation  of  the  Navier-Stokes  equations,  a  finite-difference 
scheme  was  used  to  directly  calculate  J\(ot)  and  J^a)  for  0  <  a  <  3.25.  It  is  hoped  that 
these  solutions  will  useful  as  base  line  values  for  comparison  with  other  approaches  to  these 
problems. 
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Figure  2.  Temporal  evolution  of  the  lift/drag  coefficient  at  various  values  of  speed  ratio 
(0  <  a  <  3.25 )  up  f>  t  =  24.0. 
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Figure  3.  Effect  of  speed  ratio  on  time-average  lift,  drag  and  lift/drag  coefficients  for 
0  <  a  £  3.25. 
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Figure  4.  Effect  of  speed  ratio  on  variation  of  total  lift/total  drag  force  ratio. 
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